#include<iostream>
#include<stdlib.h>
#include<iomanip>
#include<ctime>
#include <sys/time.h>
#include<arm_neon.h>
using namespace std;
const int n = 1000;//n=500,800,1000,1500,2000
float A[n][n];
void init(float A[n][n])//n*n矩阵赋值
{
    srand((unsigned)time(NULL));
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
        {
            A[i][j] = 0;
        }
        A[i][i] = 1.0;
        for (int j = i + 1; j < n; j++)
        {
            A[i][j] = rand()%10000;//模取一万保证矩阵数据不会过大
        }
    }
    for (int k = 0; k < n; k++)
    {
        for (int i = k + 1; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A[i][j] += A[k][j];
                A[i][j] = float(((int)A[i][j]) % 10000);//同理
            }
        }
    }
}
void display(float A[n][n])
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cout << setw(11) << A[i][j] << " ";
        }
        cout << endl;
    }
}
void Gauss_normal_clock(float** A)//普通高斯消元算法
{
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}
void clone(float A[n][n], float** B)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
}
void Gauss_normal_neon_one(float** A)//只优化除法部分
{
    for (int k = 0; k < n; k++)
    {
        float32x4_t vt = vdupq_n_f32(A[k][k]);
        int j;
        for (j = k + 1; j + 4 <= n; j += 4)
        {
            float32x4_t va = vld1q_f32(&A[k][j]);
            va = vdivq_f32(va, vt);
            vst1q_f32(&A[k][j], va);
        }
        for (; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}
void Gauss_normal_neon_two(float** A)//只优化减法部分
{
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            float32x4_t vaik = vdupq_n_f32(A[i][k]);
            int j;
            for (j = k + 1; j + 4 <= n; j += 4)
            {
                float32x4_t vakj = vld1q_f32(&A[k][j]);
                float32x4_t vaij = vld1q_f32(&A[i][j]);
                float32x4_t vx = vmulq_f32(vakj, vaik);
                vaij = vsubq_f32(vaij, vx);
                vst1q_f32(&A[i][j], vaij);
            }
            for (; j < n; j++)
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
}
int main()
{
    init(A);
    double time1_sum = 0;
    int count1 = 0;
    for (int l = 0; l < 10; l++)//测量平凡算法时间，运行十次求平均值
    {
        double time1 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        struct  timeval   tv_begin3, tv_end3;
        gettimeofday(&tv_begin3, NULL);
        Gauss_normal_clock(B);
        gettimeofday(&tv_end3, NULL);
        time1 = (double)(tv_end3.tv_sec - tv_begin3.tv_sec) + (tv_end3.tv_usec - tv_begin3.tv_usec) / 1000000;
        time1 *= 1000;
        if (time1 != 0)count1++;
        time1_sum += time1;
    }
    cout << "普通算法总运行时间：" << time1_sum / count1 << "ms" << endl;


    double time2_sum = 0;
    int count2 = 0;
    for (int l = 0; l < 10; l++)//测量除法优化时间，运行十次求平均值
    {
        double time2 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        struct  timeval   tv_begin3, tv_end3;
        gettimeofday(&tv_begin3, NULL);
        Gauss_normal_neon_one(B);
        gettimeofday(&tv_end3, NULL);
        time2= (double)(tv_end3.tv_sec - tv_begin3.tv_sec) + (tv_end3.tv_usec - tv_begin3.tv_usec) / 1000000;
        time2 *= 1000;
        if (time2 != 0)count2++;
        time2_sum += time2;
    }
    cout << "普通算法neon除法优化总运行时间：" << time2_sum / count2 << "ms" << endl;

    double time3_sum = 0;
    int count3 = 0;
    for(int l=0;l<10;l++)//测量减法优化时间，运行十次求平均值
    {
        double time3 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        struct  timeval   tv_begin4, tv_end4;
        gettimeofday(&tv_begin4, NULL);
        Gauss_normal_neon_two(B);
        gettimeofday(&tv_end4, NULL);
        time3 = (double)(tv_end4.tv_sec - tv_begin4.tv_sec) + (tv_end4.tv_usec - tv_begin4.tv_usec) / 1000000;
        time3 *= 1000;
        if (time3 != 0)count3++;
        time3_sum += time3;
    }
    cout << "普通算法neon减法优化总运行时间：" << time3_sum / count3 << "ms" << endl;
}
